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We study the process of heat transfer through an entangled pair of two-level system, demonstrating the role 
of quantum correlations in this nonequilibrium process. While quantum correlations generally degrade with 
increasing the temperature bias, introducing spatial asymmetry leads to an intricate behavior: Connecting the 
qubits unequally to the reservoirs one finds that quantum correlations persist and increase with the temperature 
bias when the system is more weakly linked to the hot reservoir. In the reversed case, linking the system 
more strongly to the hot bath, the opposite, more natural behavior is observed, with quantum correlations being 
strongly suppressed upon increasing the temperature bias. 
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I. INTRODUCTION 

Understanding thermal energy transfer at the nanoscale has 
recently become a topic of great interest in nanotechnology 
ifT], with proposals for new devices that can actively control 
heat conduction and information storage: thermal rectifiers 
101, thermal logic operators and memory devices iHl. As 
some of these devices have already been realized JH one 
should note that these thermal elements typically have been 
analyzed within the principles of classical mechanics, oper- 
ating at room temperature. Assessing the role of quantum 
effects in the operation of such systems is of fundamental 
and practical importance, for building quantum devices fight- 
ing relaxation and decoherence processes, operating under 
nonequilibrium conditions. 

The typical setup of interest in this context includes a 
small system, with few degrees of freedom (nanobeam, linear 
molecule, spin chain), connected at its ends to two large reser- 
voirs (solids, metals), maintained at different temperatures. In 
the steady state limit a constant heat current flows through the 
system. The thermal conductance of such a nanoscale junc- 
tion has been primarily simulated using either of the follow- 
ing three approaches: (i) The Landauer's formula J?! com- 
bined with first-principles calculations of the Hamiltonian 
force constants [8J, (ii) the non-equilibrium Green's function 
method lU i fioll . or (iii) classical molecular dynamics simu- 
lations fir]. The Boltzmann-Peierls phonon transport theory 
1 12, 13], mixed classical-quantum simulations ] 14], and quan- 
tum master equation methods ITsMitIi are other methods de- 
veloped for predicting the conductance properties of different 
objects. For systems with few degrees of freedom the latter 
method is of particular interest, as the kinetic equations of 
motion can be derived, under some approximations, from the 
fundamental quantum equations of motion, as explained be- 
low. For simple model systems these equations can be analyt- 
ically solved, providing insight on the microscopic dynamics 
IH]. 

Entanglement and quantum discord fl^ are quantum cor- 
relations with no classical counterpart which can be used as 
tools for identifying and distinguishing the quantum aspects 
in the thermal transport process from the classical ones [|l9l - 



I21I] . It is our objective here to consider a simple (yet involv- 
ing many-body interactions) thermal conducting junction, to 
study its transport characteristics within the quantum master 
equation method, and to evaluate the role of quantum effects 
in the energy transport process through the entanglement and 
discord measures. In particular, as the steady state concur- 
rence in quantum open systems has been akeady analyzed 
lfT9l - l210 . it is of interest to explore its relation to the steady 
state discord measure, quantifying non-classical correlations 
beyond entanglement. 

The particular system examined here includes a pair of a 
two-level system under a magnetic field, placed in between 
two thermal reservoirs. We use this model as a case study for 
illustrating the intricate role of the nonequilibrium condition 
on quantum correlations in the system. We analytically calcu- 
late the amount of entanglement and discord in this model at 
different baths temperatures, demonstrating that even at rel- 
atively high temperatures quantum correlations play a role 
in the energy transfer process. Furthermore, by introducing 
asymmetry, we show that, counter-intuitively, quantum corre- 
lations may be enhanced upon increasing the temperature bias 
across the system, depending on the bias polarity. 



II. MODEL 

Our model includes a quantum system coupled to two dif- 
ferent thermal reservoirs. The system incorporates two inter- 
acting qubits in a magnetic field subjected to the XY inter- 
action. These qubits, 1 and 2, are separately coupled to in- 
dependent reservoks Hl and Hji, respectively, maintained in 
thermal equilibrium at temperatures T^, v = L, i?. The total 
Hamiltonian is given by 

- + ^ + 14 + Vr, (1) 

where the two-qubit Hamiltonian is 

Hs = + + |(a^af + a^af). (2) 
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Here cr* (i = x, y, z) are the Pauli matrices. The system 
Hamiltonian can be diagonalized to produce the "system di- 
agonal" basis with four states n — 1,2, 3, 4, 
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|2) = 

|3) = 

|4> = 



-i=(|;t>-|n», E,^-K 
\U), E2 = ~e 

1 



V2 



(i;t> + in», Ei 



(3) 



The states are ordered assuming that the inter-spin interaction 
is strong, k > e > 0. As the qubits are identical, the sym- 
metric conditions UJ43 = UJ21 and W42 =^31 follow, with 
w„m — En — E„i. We assume that system-bath interactions 
are separable, and use the following bipartite form, 

K = B,S,, (4) 

with as a bath operator and Si, a system operator Specif- 
ically, we use Sl = o-f and Sr = erf . In the basis of Eq. ^ 
these operators translate to 

= -^(-|l)(2| + |l)(3| + |2)(4| + |3)(4| + /i.c.) 

= -^(|l)(2|-|l)(3| + |2)(4| + |3)(4| + /i.c.). (5) 

We next calculate the heat current and the quantum entangle- 
ment and discord. Note that we have not yet specified a model 
for the reservoirs, as the calculations can be formally done for 
different bath realizations (llTfl • 



III. STEADY STATE DYNAMICS 

Under nonequilibrium conditions, Tl 7^ Tr, in the long 
time limit, we present next the following quantities: (i) the 
system's population, (ii) the steady state heat current, and (iii) 
the resulting quantum correlations. 



A. Levels' population 

We follow standard weak coupling schemes [22] and use 
the Bom-Markov approximation. Beginning with the Liou- 
ville equation, the method first involves the assumption of 
weak system-bath interactions. Furthermore, we apply the 
Markovian limit, assuming that the reservoirs' characteristic 
timescales are shorter than the subsystem relaxation time. Un- 
der these approximations a master equation for the states pop- 
ulation Pn {n — 1,2,3,4) can be readily obtained, 

Pn{t) = l'5'mnP^r?x(^)fcm-!-n 

- Pn{t)J2\S';:,n\'K^rn- (6) 



Details about this derivation (for the two-bath case) are given 
in Refs. 1 16] and 1 17]. Here S" Si^nW^) ("I- The Fermi 
golden rule transition rates are given by 



n— »m 



dre" 



(7) 



where Bv{t) = e'^^"'^ B,je~^^"'^ are interaction picture 
operators. The thermal average is given by (0)t„ = 
Tr[e~^"'/^''0]/Tr[e~'^''/^"']. Note that the rates are evalu- 
ated at a specific subsystem frequency. For example, in Eq. 
O the relevant energy scale is tOnm = E^ — E^- Solving Eq. 
(|6]l in steady state we obtain the population 



Pi 



P2 



P^ 



Pa 



W12W13 




h W21) {Wu 


+ W31) 




W21W13 




{W12- 


h W21) (VF13 


+ W31) 




W12W31 




{W12- 


h W21) (VF13 


+ W31) 




W21W31 




{W12- 


h W21) {Wi3 


+ W31) 



(8) 



where we have introduced the short notation W,„„ = k. 



^n^m- Since the qubits are of equal energy and 

\gR 1 2 



L I 



we have also utilized the fact that Wi^ — W24 and 
W34 = 14^12 in deriving ([8]l. 

The rate constants depend on the particular choice of the 
system-bath interaction operator and the bath Hamiltonian. 
For example, assuming the reservoirs include a collection of 
harmonic modes and that the bath operator coupled to the sys- 
tem is a displacement operator. 



B,, 



i 



(9) 



the relaxation (excitation) rate with m > n (m < n)) reduces 
to 



rB,i.(a;m„)[na(w™„) + 1], 



(10) 



using the definition (Q- Here Tb,u{^) = ^Ij^i^ ~ 

uij) and nQ{uj) — [e"/-^" — 1]"-"^ is the Bose-Einstein distri- 
bution. Another physical setup is the spin reservoir, including 
a collection of P noninteracting spins. 



p=i 



B, 



p 



(11) 



Each spin is described by the two eigenstates (z = 0,1) \i)p 
and eigenenergies ep{i). In this case the relaxation rate be- 
comes 



^ S,v{,^mn)n''g{—LOmn) 



(12) 



with the spin occupation factor n.^(a;) = [e'^''^" + 1] ^ 

and the effective spin-bath-system coupling r5.i,(w) — 
2 

5{uj + ep(0) — ep(l)). For details see 



Ref. ill. 
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B. Heat current 

Formally, the expectation value of the current, calculated, 
e.g., at the left contact, is given by Jl ~ ^Tr{[HL — 
Hs: Vl]p), P is the total density matrix ll23ll . In steady state, 
the expectation value of the interaction is zero, Tr(^^p) — 
iTi([HL + Hs, Vl]p), and this expression reduces to Jl = 
iTr([VL, Hs]p)- Using the system-diagonal representation, it 
is strai ghtf orward to show that the steady state current be- 
comes ilTIl 



(13) 




where w„j„ — Em — En. Under the Born-Markov approxima- 
tion used above to resolve the population dynamics, a second- 
order expression for the steady-state current can be further ob- 
tained 01 



Jl = 



\S,. 



L I 

mn 



(14) 



FIG. 1: Concurrence (full) and discord (dashed) in an equilibrium 
system, Ta = Tl ~ Tr, The inset displays the levels population, top 
(Pi) to bottom (P4). The physical parameters are e = 0.2, k = 1. 
Bosonic reservoirs were adopted. 



D. Discord 



We apply this expression on the two-qubit model. A some- 
what tedious calculation shows that the current is given by a 
sum of two terms 
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(15) 



C. Concurrence 

The entanglement of formation is a monotonically increas- 
ing function of the Wootters' concurrence [24J. We calcu- 
late next the nonequilibrium concurrence in the two-qubit 
model. The states |1) and |4) are entangled, and the na- 
ture of the entanglement may remain in the final steady state, 
which is described by the diagonal reduced density matrix 
Pd = diag{Pi, P2, P3, P4), in the eigenbasis of the system 
Hamiltonian Hs- On the other hand, in the uncoupled basis 
Jit) , Itl) and Itt), the reduced density matrix is given 
by a nondiagonal form. 
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(16) 



Using the general form given in, e.g., [25], the concurrence 
can be expressed in the uncoupled basis by. 



C{Tl, Tr) = max(2P„,ax - Pi - P4 - 2^/^, 0) (17) 



where P„iax — max(Pi, P4, v^^)- This function in gen- 
eral depends on the temperatures of both thermal baths. 



Quantum discord ifTsll . quantifying nonlocal correlations, is 
given by the difference between the quantum mutual informa- 
tion I{p) and the classical correlation C(p), 



Q{p)^I{p)-C{p). 



(18) 



The analytical expressions for the classical correlation and the 
quantum discord were obtained for a class of X states in Ref. 
i25ll . Adjusting these expression to our model, the quantum 
mutual information can be written in terms of the steady state 
populations ^ as 

Ap) = 

2 - log2[(l - P2 + P^Y-P-+P-{1 +P2- P3)l+^^-^3] 
+Pl log2 Pi + P2 l0g2 P2 + P3 log2 P3 + P4 l0g2 P4, 

(19) 

whereas the classical correlation is given by 



C(p) = 
1 
2 

min{5i, 5*2} 



1 - o log2[(l -P2 + P3)'-^^+^^(l + P2- P^f 



(20) 



Here 
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and 



52 = l--l0g2[(l-X)^-^(l + /f)^ + ^]. 



(22) 
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FIG. 2: (a) Nonequilibrium thermal correlations for a pair of qubits 
coupled to boson baths. Concurrence (full) and discord (dashed) for 
Tl = 1.5, Tr is modified, (b) Same for Tl = 0.2. The two-qubit 
parameters are e = 0.2, k = \ and we set Vb,v = 1. 
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FIG. 3: (a) Nonequilibrium thermal correlations for a pair of qubits 
coupled to spin baths. Concurrence (full) and discord (dashed) for 
Tl = 1.5, Tr is modified, (b) Same for Tl = 0.2. The two-qubit 
parameters are e = 1 and k = 0.2, and we set Vs,v = 1. 



The coefficient K is defined 



as 



K 



IV. EXAMPLES 

In what follows we analyze the quantum correlations, con- 
currence and discord, for various nonequilibrium conditions, 
showing that quantum discord survives at relatively high tem- 
peratures where concurrence is zero. We also study an asym- 
metric scenario demonstrating that quantum correlations may 
exist even at large temperature biases. We will typically use 
the following parameters: spin energy e — 0.2 and large spin- 
spin interaction k = 1. We will also assume that the relaxation 
rates do not depend on energy, thus treat them as constants, 
Tb^u and Ts.iy The factor l/\/2 in Eq. (|5]l is absorbed into 



the definition of the rates T. 

We begin our analysis with an equilibrium situation, Ta = 
Tl — Tr, and compare the concurrence and the discord mea- 
sures at different temperatures. Fig. [T]shows that at very low 
temperatures both measures yield the same result. At higher 
temperatures (yet Ta < k) concurrence is slightly larger than 
discord. At even higher temperatures, Ta > k, concurrence 
suddenly diminishes [26] while quantum discord is still finite, 
slowly decreasing to zero. This is in accord with the fact that 
discord quantifies non-classical correlations beyond entangle- 
ment. Generally, both quantum correlations decay with tem- 
perature due to thermal relaxation effects. As can be inferred 
from the inset, when the ground state population falls below 
~ 1/2 the concurrence dies. We also found that these ob- 
servations were not sensitive to the reservoirs propeities, and 
similar trends were obtained using either boson or spin baths. 

Fig. ID displays the nonequilibrium thermal correlations, 
keeping Tl fixed and changing T/j. In the classical limit, 
at high temperatures (top panel) the discord overcomes the 
concurrence, even when Tn is low. In contrast, at low tem- 
peratures (bottom panel) the opposite trend is observed till a 
crossover value beyond which the concurrence dies ||26| yet 
the discord is finite. While in Fig. |2]bosonic reservoirs have 
been adopted. Fig. [3]displays the nonequilibrium thermal cor- 
relations using spin baths. The main difference noted is that 
when Tl is quite high and is low [see Fig. Oa)], the spin- 
bath case shows that concuiTence is higher than discord, while 
the opposite behavior is observed in Fig. I2a). This obser- 
vation is in accord with the notion that spin reservoirs are 
"more quantum" than harmonic baths in the sense that har- 
monic modes can be represented by two-level systems (spins) 
at low enough temperatures. 

Next, we explore the role of spatial asymmetry on the sur- 
vival of quantum correlations in (temperature) driven systems. 
Thermal rectification, an asymmetry of the heat current for 
forward and reversed temperature gradients, has been exten- 
sively analyzed in the last decade iHIsI ITtIi . In a desirable rec- 
tifier the system behaves as an excellent heat conductor in one 
direction of the temperature bias, while for the opposite direc- 
tion it effectively acts as an insulator It is agreed that junc- 
tions incorporating anharmonic interactions with some sort of 
spatial asymmetry should demonstrate this effect. The two- 
qubit model, prepared with some asymmetry, e.g., assuming 
that the qubits asymmetrically couple to, e.g., bosonic reser- 
voirs, Tb,l 7^ ^B,R, is expected to behave as a thermal recti- 
fier. 

Figure |4] indeed shows the emergence of the thermal rec- 
tifying effect upon turning on the asymmetry. The current 
in Fig. SJb) is symmetric since Tb,l = ^b.r- In con- 
trast, in Fig. Ud) the heat current is larger (in magnitude) 
when the system is more strongly linked to the cold reservoir 
(Tl < Tr and Tb.l > ^b,fi) Ii2|. Here the averaged temper- 
ature is Ta = 1 with Tl ^ Ta + AT and Tr ^ Ta - AT. 
While the effect of thermal rectification is well understood, 
here we demonstrate that the transition between the fairly con- 
ducting phase (AT < 0) and the poorly conducting phase 
(AT > 0) corresponds to a turnover in the transport mech- 
anism. In the symmetric setup discord and concurrence are 
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FIG. 4: Concurrence (full) and discord (dashed) in symmetric (a) 
and asymmetric systems (c), Ta = 1 is fixed and the temperatures at 
the two ends are modulated, Tl = Ta + AT, = Ta - AT. In 
the symmetric case Tb,l ~ ^b.r = 1. Asymmetiy is introduced 
by taking Tb.l/Tb.r = 20 with Tb.l = 1 Panels (b) and (d) 
further display the current for the symmetric and asymmetric cases, 
respectively. In all panels e = 0.2, k = 1. 



symmetric functions, see Fig. g^a), C(AT) = C(-AT) and 
Q(AT) = Q(-AT). Fig. I^c) shows a different behavior in 
the presence of asymmetry: When the bias is negative quan- 
tum correlations persist. However, for the opposite polarity 
the concurrence is zero and discord is diminishing. 

We can reason this behavior by noting that when the sys- 
tem is more strongly attached to the cold bath, Fs ^ > Fs j^ 
and Tr > Tl, the ground state population is larger than that 
expected in the opposite case. Since the ground state is an en- 
tangled (singlet) state, |1) = 75 (lit) - Iti). its large popu- 
lation reflects energy transmission process assisted by a quan- 
tum correlated system state. 

By further increasing k one can extend the range over which 
quantum correlations survive, see Fig. |5] For the asymmetric 
setup again we note that when T^ < Tji the current is large, 
in comparison to the opposite polarity, and that entanglement 
measures are close to unity. For the reversed case, Tl > Tn, 
the current magnitude is lowered, and quantum correlations 
are being suppressed. In particular, for Tl—T/j = —2, discord 
is large, Q ^ 1, while at Tl — Tr = 2 it is reduced by an order 
of magnitude, Q ^ 0.1. Another interesting observation is 
that both discord and concurrence display a non-monotonic 
behavior at large negative bias, reflected in a small maxima 
around AT 1. 7 1 19-21]. 

To conclude this section, by switching the sign of the tem- 
perature bias one can control the magnitude of the heat current 
in asymmetric spin chains and the underlying transport mech- 
anism, as reflected by the survival or suppression of quantum 
correlations in the system. 
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FIG. 5: Concurrence (full) and discord (dashed) in symmetric (a) 
and asymmetric systems (c). Panels (b) and (d) display the current for 
the symmetric and asymmetric situations, respectively. Parameters 
are the same as in Fig. |4l but the inter-spin interaction has been 
increased to k = 2. 



V. SUMMARY 

We detailed here a simple model of a many-body open 
quantum system which could be analytically solved, useful 
for analyzing the role of the temperature gradient on quantum 
correlations in a conducting nanojunction. Our calculations 
manifest that for symmetric systems under large temperature 
gradients and at high temperatures quantum discord can be 
maintained, slowly decaying with AT, whereas concurrence 
typically dies. In the presence of asymmetry we found that 
quantum correlations can survive in one direction of the tem- 
perature gradient, while they diminish when reversing the bias 
direction. 

The present analysis could be generalized for describ- 
ing transport and quantum correlations in longer-linear spin 
chains. One could also treat other systems, e.g., large spins, or 
adopt unequal reservoirs at the two ends, for example, assum- 
ing the system is coupled to both a solid and a metal 1 17]. By 
complementing transport studies with the calculation of quan- 
tum correlations one can estimate and corroborate the role of 
quantum effects in the transport process. This might be use- 
ful for building quantum devices operating in noisy-thermal 
environments under nonequilibrium conditions iIztIi . 
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